#include "numuCCMultiPiAnalysis.hxx"
#include "Parameters.hxx"
#include "CategoriesUtils.hxx"
#include "BasicUtils.hxx"


const unsigned int NMAXTPCSEC = 20;
const unsigned int NMAXFGDSEC = 20;
const unsigned int NMAXTPC1 = 20;

//********************************************************************
numuCCMultiPiAnalysis::numuCCMultiPiAnalysis(AnalysisAlgorithm* ana) : baseAnalysis(ana){
//********************************************************************

  // Minimum accum level to save event into the output tree
  SetMinAccumCutLevelToSave(ND::params().GetParameterI("numuCCMultiPiAnalysis.MinAccumLevelToSave"));

  // Create a numuCCAnalysis passing this analysis to the constructor. In that way the same managers are used
  _numuCCAnalysis = new numuCCAnalysis(this);

  // Use the numuCCAnalysis (in practice that means that the sabe box and event will be used for the numuCCAnalysis as for this analysis)
  // By default configurations, systematics, corrections and tree variables of the used algorithm will be added
  UseAnalysis(_numuCCAnalysis);

  // Add the package version
  ND::versioning().AddPackage("numuCCMultiPiAnalysis", anaUtils::GetSoftwareVersionFromPath((std::string)getenv("NUMUCCMULTIPIANALYSISROOT")));
}

//********************************************************************
void numuCCMultiPiAnalysis::DefineSelections(){
//********************************************************************

  // Add the selection to the Selection Manager
  sel().AddSelection("kTrackerNumuCCMultiPi",    "numuCC multi-pionselection",     new numuCCMultiPiSelection());
}

//********************************************************************
void numuCCMultiPiAnalysis::DefineMicroTrees(bool addBase){
//********************************************************************

  // -------- Add variables to the analysis tree ----------------------
  
  // Variables from the numuCCAnalysis analysis (including the ones in baseAnalysis by default, otherwise addBase should be false
  if (addBase) _numuCCAnalysis->DefineMicroTrees(addBase);

  // --- Add here more variables for the CCMultiPi analysis -----

  AddVarI(output(),MuonVertexId, "");
  AddVarI(output(),MuonIndex,"");

  // --- Pion candidates 
  AddVarVF(output(),NegPionMom,"",       NNegPion);
  AddVarVF(output(),NegPionPidLik,"",    NNegPion);
  AddVarVF(output(),NegPionElPidLik,"",  NNegPion);
  AddVarVI(output(),NegPionTId,"",       NNegPion);
  AddVarVI(output(),NegPionIndex,"",     NNegPion);
  AddVarVI(output(),NegPionParentTId,"", NNegPion);
  AddVarVI(output(),NegPionGParentTId,"",NNegPion);
  AddVarVI(output(),NegPionVId,"",       NNegPion);
  AddVarMF(output(),NegPionDir,"",       NNegPion,-NMAXTPCSEC,3);
  AddVarMF(output(),NegPionPosStart,"",  NNegPion,-NMAXTPCSEC,3);
  AddVarMF(output(),NegPionPosEnd,"",    NNegPion,-NMAXTPCSEC,3);
  
  AddVarVF(output(),PosPionMom,"",       NPosPion);
  AddVarVF(output(),PosPionPidLik,"",    NPosPion);
  AddVarVF(output(),PosPionElPidLik,"",  NPosPion);
  AddVarVI(output(),PosPionTId,"",       NPosPion);
  AddVarVI(output(),PosPionIndex,"",     NPosPion);
  AddVarVI(output(),PosPionParentTId,"", NPosPion);
  AddVarVI(output(),PosPionGParentTId,"",NPosPion);
  AddVarVI(output(),PosPionVId,"",       NPosPion);
  AddVarMF(output(),PosPionDir,"",       NPosPion,-NMAXTPCSEC,3);
  AddVarMF(output(),PosPionPosStart,"",  NPosPion,-NMAXTPCSEC,3);
  AddVarMF(output(),PosPionPosEnd,"",    NPosPion,-NMAXTPCSEC,3);
  
  // --- Pi0 candidates
  AddVarVF(output(),Pi0ElMom,"",       NPi0El);
  AddVarVF(output(),Pi0ElPull,"",      NPi0El);
  AddVarVI(output(),Pi0ElTId,"",       NPi0El);
  AddVarVI(output(),Pi0ElIndex,"",     NPi0El);
  AddVarVI(output(),Pi0ElParentTId,"", NPi0El);
  AddVarVI(output(),Pi0ElGParentTId,"",NPi0El);
  AddVarVI(output(),Pi0ElVId,"",       NPi0El);
  AddVarMF(output(),Pi0ElDir,"",       NPi0El,-NMAXTPCSEC,3);
  AddVarMF(output(),Pi0ElPosStart,"",  NPi0El,-NMAXTPCSEC,3);
  AddVarMF(output(),Pi0ElPosEnd,"",    NPi0El,-NMAXTPCSEC,3);

  AddVarVF(output(),Pi0PosMom,"",       NPi0Pos);
  AddVarVF(output(),Pi0PosPull,"",      NPi0Pos);
  AddVarVI(output(),Pi0PosTId,"",       NPi0Pos);
  AddVarVI(output(),Pi0PosIndex,"",     NPi0Pos);
  AddVarVI(output(),Pi0PosParentTId,"", NPi0Pos);
  AddVarVI(output(),Pi0PosGParentTId,"",NPi0Pos);
  AddVarVI(output(),Pi0PosVId,"",       NPi0Pos);
  AddVarMF(output(),Pi0PosDir,"",       NPi0Pos,-NMAXTPCSEC,3);
  AddVarMF(output(),Pi0PosPosStart,"",  NPi0Pos,-NMAXTPCSEC,3);
  AddVarMF(output(),Pi0PosPosEnd,"",    NPi0Pos,-NMAXTPCSEC,3);
 
  // --- ALL TPC sec tracks 
  AddVarVF(output(),TPCSecMom,"",       NTPCSec);
  AddVarVF(output(),TPCSecMuPidLik,"",  NTPCSec);
  AddVarVF(output(),TPCSecPiPidLik,"",  NTPCSec);
  AddVarVF(output(),TPCSecPrPidLik,"",  NTPCSec);
  AddVarVF(output(),TPCSecElPidLik,"",  NTPCSec);
  AddVarVI(output(),TPCSecDetectors,"", NTPCSec);
  AddVarVI(output(),TPCSecQ,"",         NTPCSec);
  AddVarVI(output(),TPCSecTId,"",       NTPCSec);
  AddVarVI(output(),TPCSecIndex,"",     NTPCSec);
  AddVarVI(output(),TPCSecParentTId,"", NTPCSec);
  AddVarVI(output(),TPCSecGParentTId,"",NTPCSec);
  AddVarMF(output(),TPCSecDir,"",       NTPCSec,-NMAXTPCSEC,3);
  AddVarMF(output(),TPCSecPosStart,"",  NTPCSec,-NMAXTPCSEC,3);
  AddVarMF(output(),TPCSecPosEnd,"",    NTPCSec,-NMAXTPCSEC,3);
  
  // ---- IsoFGD Pion tracks -----
  AddVarVF(output(),FGDPiLength,"",           NFGDPi); 
  AddVarVF(output(),FGDPiPiPull,"",           NFGDPi); 
  AddVarVF(output(),FGDPiAvgTime,"",          NFGDPi); 
  AddVarVF(output(),FGDPiDistance,"",         NFGDPi);
  AddVarVF(output(),FGDPiMuonangle,"",        NFGDPi);  
  AddVarVF(output(),FGDPiCosTheta,"",         NFGDPi);
  AddVarVI(output(),FGDPiTId,"",              NFGDPi);
  AddVarVI(output(),FGDPiIndex,"",            NFGDPi);
  AddVarVI(output(),FGDPiParentTId,"",        NFGDPi);
  AddVarVI(output(),FGDPiGParentTId,"",       NFGDPi);
  AddVarVI(output(),FGDPiVId,"",              NFGDPi);
  AddVarMF(output(),FGDPiDir,"",              NFGDPi,-NMAXTPCSEC,3);
  AddVarMF(output(),FGDPiPosStart,"",         NFGDPi,-NMAXTPCSEC,3);
  AddVarMF(output(),FGDPiPosEnd,"",           NFGDPi,-NMAXTPCSEC,3);
  AddVarVF(output(),FGDPiTimeIni,"",          NFGDPi); 
  AddVarVF(output(),FGDPiTimeEnd,"",          NFGDPi); 
  AddVarVF(output(),FGDPiDeltaTimeIniSelMu,"",NFGDPi);  
  AddVarVF(output(),FGDPiDeltaTimeEndSelMu,"",NFGDPi); 

  // ---- IsoFGD electron and positron tracks -----
  AddVarVF(output(),FGDElLength,"",           NFGDEl); 
  AddVarVF(output(),FGDElPiPull,"",           NFGDEl); 
  AddVarVF(output(),FGDElAvgTime,"",          NFGDEl); 
  AddVarVI(output(),FGDElTId,"",              NFGDEl); 
  AddVarVI(output(),FGDElIndex,"",            NFGDEl); 
  AddVarVI(output(),FGDElParentTId,"",        NFGDEl); 
  AddVarVI(output(),FGDElGParentTId,"",       NFGDEl); 
  AddVarVI(output(),FGDElVId,"",              NFGDEl);
  AddVarVF(output(),FGDElTimeIni,"",          NFGDEl);  
  AddVarVF(output(),FGDElTimeEnd,"",          NFGDEl); 
  AddVarVF(output(),FGDElDeltaTimeIniSelMu,"",NFGDEl);  
  AddVarVF(output(),FGDElDeltaTimeEndSelMu,"",NFGDEl); 
  AddVarVF(output(),FGDElDistance,"",         NFGDEl);
  AddVarVF(output(),FGDElMuonangle,"",        NFGDEl); 
  AddVarVF(output(),FGDElCosTheta,"",         NFGDEl);
  AddVarMF(output(),FGDElDir,"",              NFGDEl,-NMAXFGDSEC,3);
  AddVarMF(output(),FGDElPosStart,"",         NFGDEl,-NMAXFGDSEC,3);
  AddVarMF(output(),FGDElPosEnd,"",           NFGDEl,-NMAXFGDSEC,3);


  // ---- IsoFGD tracks -----
  AddVarVF(output(),FGDSecLength,"",           NFGDSec); 
  AddVarVF(output(),FGDSecPiPull,"",           NFGDSec); 
  AddVarVF(output(),FGDSecAvgTime,"",          NFGDSec); 
  AddVarVI(output(),FGDSecContained,"",        NFGDSec); 
  AddVarVI(output(),FGDSecTId,"",              NFGDSec); 
  AddVarVI(output(),FGDSecIndex,"",            NFGDSec); 
  AddVarVI(output(),FGDSecParentTId,"",        NFGDSec); 
  AddVarVI(output(),FGDSecGParentTId,"",       NFGDSec);
  AddVarVF(output(),FGDSecTimeIni,"",          NFGDSec); 
  AddVarVF(output(),FGDSecTimeEnd,"",          NFGDSec); 
  AddVarVF(output(),FGDSecDeltaTimeIniSelMu,"",NFGDSec);  
  AddVarVF(output(),FGDSecDeltaTimeEndSelMu,"",NFGDSec); 
  AddVarVF(output(),FGDSecCosTheta,"",         NFGDSec);
  AddVarMF(output(),FGDSecDir,"",              NFGDSec,-NMAXFGDSEC,3);
  AddVarMF(output(),FGDSecPosStart,"",         NFGDSec,-NMAXFGDSEC,3);
  AddVarMF(output(),FGDSecPosEnd,"",           NFGDSec,-NMAXFGDSEC,3);
 
  
  // ---- Michel Eletrons ----
  AddVarVI(output(),MENHits,"",    NME); 
  AddVarVF(output(),MERawCharge,"",NME); 
  AddVarVF(output(),MEMinTime,"",  NME); 
  AddVarVF(output(),MEMaxTime,"",  NME); 
  
  // ---- TPC1 Tracks ----
  AddVarVF(output(),TPC1TrackMom,"",       NTPC1Track); 
  AddVarVF(output(),TPC1TrackCosTheta,"",  NTPC1Track);
  AddVarVF(output(),TPC1TrackPhi,"",       NTPC1Track);
  AddVarVI(output(),TPC1TrackVId,"",       NTPC1Track); 
  AddVarVI(output(),TPC1TrackTId,"",       NTPC1Track); 
  AddVarVI(output(),TPC1TrackParentTId,"", NTPC1Track); 
  AddVarVI(output(),TPC1TrackGParentTId,"",NTPC1Track); 
  AddVarMF(output(),TPC1TrackPosStart,"",  NTPC1Track,-NMAXTPC1,3); 
  AddVarMF(output(),TPC1TrackPosEnd,"",    NTPC1Track,-NMAXTPC1,3); 
 
}

//********************************************************************
void numuCCMultiPiAnalysis::FillMicroTrees(bool addBase){
//********************************************************************

  // Variables from the numuCCAnalysis analysis (including the ones in baseAnalysis by default, otherwise addBase should be false
  if (addBase) _numuCCAnalysis->FillMicroTrees(addBase);
  
  if (mybox().HMNtrack  ) {
    if( mybox().HMNtrack->TrueTrack ){
      output().FillVar(MuonIndex,mybox().HMNtrack->TrueTrack->ID);
      if( mybox().HMNtrack->TrueTrack->TrueVertex ) 
	output().FillVar(MuonVertexId,static_cast<AnaTrueVertex*>(mybox().HMNtrack->TrueTrack->TrueVertex)->RooVtxIndex);
    }
  }

  // --- Fill here more variables for the CCMultiPi analysis -----
  // Selected negative pions 
  for( Int_t i = 0; i < mybox().nNegativePionTPCtracks; i++ ){
    AnaTrackB *track = mybox().NegativePionTPCtracks[i];
    Float_t ellklh = (anaUtils::GetPIDLikelihood(*track,0)+anaUtils::GetPIDLikelihood(*track,3))/(1.-anaUtils::GetPIDLikelihood(*track,2)); 
    output().FillVectorVar(NegPionMom,      track->Momentum);
    output().FillVectorVar(NegPionPidLik,   anaUtils::GetPIDLikelihood(*track,3));    
    output().FillVectorVar(NegPionPidLik,   ellklh);
    if( track->TrueTrack ){
      output().FillVectorVar(NegPionTId,  track->TrueTrack->PDG);
      output().FillVectorVar(NegPionIndex,track->TrueTrack->ID); 
      output().FillVectorVar(NegPionParentTId, track->TrueTrack->ParentPDG);
      output().FillVectorVar(NegPionGParentTId,track->TrueTrack->GParentPDG);
      if( track->TrueTrack->TrueVertex )
	output().FillVectorVar(NegPionVId,static_cast<AnaTrueVertex*>(track->TrueTrack->TrueVertex)->RooVtxIndex);
    }
    output().FillMatrixVarFromArray(NegPionPosStart, track->PositionStart, 3);
    output().FillMatrixVarFromArray(NegPionPosEnd,   track->PositionEnd,3);
    output().FillMatrixVarFromArray(NegPionDir,      track->DirectionStart,3);
    output().IncrementCounter(NNegPion);
  }

  // Selected positive pions 
  for( Int_t i = 0; i <mybox().nPositivePionTPCtracks; i++ ){
    AnaTrackB *track = mybox().PositivePionTPCtracks[i];
    Float_t ellklh = (anaUtils::GetPIDLikelihood(*track,0)+anaUtils::GetPIDLikelihood(*track,3))/(1.-anaUtils::GetPIDLikelihood(*track,2)); 
    output().FillVectorVar(PosPionMom,      track->Momentum);
    output().FillVectorVar(PosPionPidLik,   anaUtils::GetPIDLikelihood(*track,3));    
    output().FillVectorVar(PosPionPidLik,   ellklh);

    if( track->TrueTrack ){
      output().FillVectorVar(PosPionTId,  track->TrueTrack->PDG);
      output().FillVectorVar(PosPionIndex,track->TrueTrack->ID); 
      output().FillVectorVar(PosPionParentTId, track->TrueTrack->ParentPDG);
      output().FillVectorVar(PosPionGParentTId,track->TrueTrack->GParentPDG);
      if( track->TrueTrack->TrueVertex )
	output().FillVectorVar(PosPionVId,static_cast<AnaTrueVertex*>(track->TrueTrack->TrueVertex)->RooVtxIndex);
    }

    output().FillMatrixVarFromArray(PosPionPosStart, track->PositionStart, 3);
    output().FillMatrixVarFromArray(PosPionPosEnd,   track->PositionEnd,3);
    output().FillMatrixVarFromArray(PosPionDir,      track->DirectionStart,3);
    output().IncrementCounter(NPosPion);
  }
    
  // Selected pi0 electron tracks 
  for( Int_t i = 0; i < mybox().nElPi0TPCtracks; i++ ) {
    AnaTrackB *track = mybox().ElPi0TPCtracks[i];
    output().FillVectorVar(Pi0ElMom, track->Momentum);
    output().FillVectorVar(Pi0ElPull, anaUtils::GetPIDLikelihood(*track,1)); 
   if( track->TrueTrack ){
      output().FillVectorVar(Pi0ElTId,       track->TrueTrack->PDG);
      output().FillVectorVar(Pi0ElIndex,     track->TrueTrack->ID);
      output().FillVectorVar(Pi0ElParentTId, track->TrueTrack->ParentPDG);
      output().FillVectorVar(Pi0ElGParentTId,track->TrueTrack->GParentPDG);
      if( track->TrueTrack->TrueVertex  )
	output().FillVectorVar(Pi0ElVId,static_cast<AnaTrueVertex*>(track->TrueTrack->TrueVertex)->RooVtxIndex);

    } 

    output().FillMatrixVarFromArray(Pi0ElPosStart, track->PositionStart, 3);
    output().FillMatrixVarFromArray(Pi0ElPosEnd,   track->PositionEnd,   3);
    output().FillMatrixVarFromArray(Pi0ElDir,      track->DirectionStart,3);
    output().IncrementCounter(NPi0El);

  }

  // Selected pi0 positron tracks 
  for( Int_t i = 0; i < mybox().nPosPi0TPCtracks; i++ ) {
    AnaTrackB *track = mybox().PosPi0TPCtracks[i];
    output().FillVectorVar(Pi0PosMom, track->Momentum);
    output().FillVectorVar(Pi0PosPull, anaUtils::GetPIDLikelihood(*track,1)); 
   if( track->TrueTrack ){
      output().FillVectorVar(Pi0PosTId,       track->TrueTrack->PDG);
      output().FillVectorVar(Pi0PosIndex,     track->TrueTrack->ID);
      output().FillVectorVar(Pi0PosParentTId, track->TrueTrack->ParentPDG);
      output().FillVectorVar(Pi0PosGParentTId,track->TrueTrack->GParentPDG);
      if( track->TrueTrack->TrueVertex  )
	output().FillVectorVar(Pi0PosVId,static_cast<AnaTrueVertex*>(track->TrueTrack->TrueVertex)->RooVtxIndex);

    } 

    output().FillMatrixVarFromArray(Pi0PosPosStart, track->PositionStart, 3);
    output().FillMatrixVarFromArray(Pi0PosPosEnd,   track->PositionEnd,   3);
    output().FillMatrixVarFromArray(Pi0PosDir,      track->DirectionStart,3);
    output().IncrementCounter(NPi0Pos);

  }
 

  // All TPC secondary tracks 
  for( Int_t i = 0; i < mybox().nPositiveTPCtracks; i++ ) {
    AnaTrackB *track = mybox().PositiveTPCtracks[i];

    if( mybox().HMNtrack == track ) continue; // This is the muon

    output().FillVectorVar(TPCSecMom,               track->Momentum);
    output().FillVectorVar(TPCSecPiPidLik,          anaUtils::GetPIDLikelihood(*track,3));
    output().FillVectorVar(TPCSecMuPidLik,          anaUtils::GetPIDLikelihood(*track,0));
    output().FillVectorVar(TPCSecElPidLik,          anaUtils::GetPIDLikelihood(*track,1));
    output().FillVectorVar(TPCSecPrPidLik,          anaUtils::GetPIDLikelihood(*track,2));
    output().FillMatrixVarFromArray(TPCSecPosStart, track->PositionStart,3);
    output().FillMatrixVarFromArray(TPCSecPosEnd,   track->PositionEnd,3);
    output().FillMatrixVarFromArray(TPCSecDir,      track->DirectionStart,3);
    output().FillVectorVar(TPCSecQ,1); 
    if( track->TrueTrack ){
      output().FillVectorVar(TPCSecTId,       track->TrueTrack->PDG);
      output().FillVectorVar(TPCSecIndex,     track->TrueTrack->ID);
      output().FillVectorVar(TPCSecParentTId, track->TrueTrack->ParentPDG);
      output().FillVectorVar(TPCSecGParentTId,track->TrueTrack->GParentPDG);
    }
    
    output().IncrementCounter(NTPCSec); 
  }


   for( Int_t i = 0; i < mybox().nNegativeTPCtracks; i++ ) {
    AnaTrackB *track = mybox().NegativeTPCtracks[i];
    if( mybox().HMNtrack == track ) continue; // This is the muon

    output().FillVectorVar(TPCSecMom,               track->Momentum);
    output().FillVectorVar(TPCSecPiPidLik,          anaUtils::GetPIDLikelihood(*track,3));
    output().FillVectorVar(TPCSecMuPidLik,          anaUtils::GetPIDLikelihood(*track,0));
    output().FillVectorVar(TPCSecElPidLik,          anaUtils::GetPIDLikelihood(*track,1));
    output().FillVectorVar(TPCSecPrPidLik,          anaUtils::GetPIDLikelihood(*track,2));
    output().FillMatrixVarFromArray(TPCSecPosStart, track->PositionStart,3);
    output().FillMatrixVarFromArray(TPCSecPosEnd,   track->PositionEnd,3);
    output().FillMatrixVarFromArray(TPCSecDir,      track->DirectionStart,3);
    output().FillVectorVar(TPCSecQ,-1); 
    if( track->TrueTrack ){
      output().FillVectorVar(TPCSecTId,       track->TrueTrack->PDG);
      output().FillVectorVar(TPCSecIndex,     track->TrueTrack->ID);
      output().FillVectorVar(TPCSecParentTId, track->TrueTrack->ParentPDG);
      output().FillVectorVar(TPCSecGParentTId,track->TrueTrack->GParentPDG);
    }
    output().IncrementCounter(NTPCSec); 
  }
  

  // isoFGD pion candidates 
  for (Int_t i=0;i < mybox().nIsoFGDPiontracks;i++){
    AnaTrackB* track =  mybox().IsoFGDPiontracks[i];
    AnaFgdTrack *fgdTrack = static_cast<AnaFgdTrack*>(track->FGDSegments[0]);
    if( !fgdTrack ) continue; 
    ///checking times
    Float_t timeinipi = fgdTrack->PositionStart[3];
    Float_t timeendpi = fgdTrack->PositionEnd[3];

    output().FillVectorVar(FGDPiPiPull,  fgdTrack->Pullpi); 
    output().FillVectorVar(FGDPiAvgTime, fgdTrack->AvgTime); 
    // Direction of the segment 
    output().FillVectorVar(FGDPiCosTheta,(Float_t)(utils::ArrayToTVector3(fgdTrack->DirectionStart).CosTheta()));
    //output().FillVectorVar(FGDPiLength,fgdTrack->Length); 
    output().FillVectorVar(FGDPiLength,  fgdTrack->X); 
    output().FillVectorVar(FGDPiTimeIni, timeinipi);  
    output().FillVectorVar(FGDPiTimeEnd, timeendpi); 
    
    output().FillMatrixVarFromArray(FGDPiDir,     fgdTrack->DirectionStart,3);
    output().FillMatrixVarFromArray(FGDPiPosStart,fgdTrack->PositionStart, 3);
    output().FillMatrixVarFromArray(FGDPiPosEnd,  fgdTrack->PositionEnd,   3);
    if( track->TrueTrack ){
      output().FillVectorVar(FGDPiTId,       track->TrueTrack->PDG);
      output().FillVectorVar(FGDPiIndex,     track->TrueTrack->ID); 
      output().FillVectorVar(FGDPiParentTId, track->TrueTrack->ParentPDG);
      output().FillVectorVar(FGDPiGParentTId,track->TrueTrack->GParentPDG);
      if( track->TrueTrack->TrueVertex )
      	output().FillVectorVar(FGDPiVId,static_cast<AnaTrueVertex*>(track->TrueTrack->TrueVertex)->RooVtxIndex);
     }
     if( mybox().HMNtrack ) { 
       Float_t deltainimuinipi = TMath::Abs(timeinipi-mybox().HMNtrack->PositionStart[3]);
       Float_t deltaendmuinipi = TMath::Abs(timeendpi-mybox().HMNtrack->PositionStart[3]); 
       output().FillVectorVar(FGDPiDeltaTimeIniSelMu,deltainimuinipi);  
       output().FillVectorVar(FGDPiDeltaTimeEndSelMu,deltaendmuinipi); 
       // Minimum distance to the muon vertex 
       Float_t dist1 = (utils::ArrayToTVector3(fgdTrack->PositionStart) - utils::ArrayToTVector3(mybox().HMNtrack->PositionStart)).Mag();
       Float_t dist2 = (utils::ArrayToTVector3(fgdTrack->PositionEnd)   - utils::ArrayToTVector3(mybox().HMNtrack->PositionStart)).Mag();
       output().FillVectorVar(FGDPiDistance,TMath::Min(dist1,dist2)); 
     }
     output().IncrementCounter(NFGDPi);
  }

  // isoFGD electron and positron candidates 
  for (Int_t i=0;i < mybox().nIsoFGDElPi0tracks;i++){
    AnaTrackB* track =  mybox().IsoFGDElPi0tracks[i];
    AnaFgdTrack *fgdTrack = static_cast<AnaFgdTrack*>(track->FGDSegments[0]);
    if( !fgdTrack ) continue; 
    Float_t timeiniEl = fgdTrack->PositionStart[3];
    Float_t timeendEl = fgdTrack->PositionEnd[3];
    output().FillVectorVar(FGDElCosTheta,(Float_t)(utils::ArrayToTVector3(fgdTrack->DirectionStart).CosTheta()));
    output().FillVectorVar(FGDElPiPull,fgdTrack->Pullpi); 
    output().FillVectorVar(FGDElAvgTime,fgdTrack->AvgTime); 
    //output().FillVectorVar(FGDElLength,fgdTrack->Length); 
    output().FillVectorVar(FGDElLength,fgdTrack->X); 
    output().FillVectorVar(FGDElTimeIni,timeiniEl);  
    output().FillVectorVar(FGDElTimeEnd,timeendEl); 
    //  Minimum distance to the muon vertex 
    if( mybox().HMNtrack ) {
      Float_t deltainimuiniEl = TMath::Abs(timeiniEl-mybox().HMNtrack->PositionStart[3]);
      Float_t deltaendmuiniEl = TMath::Abs(timeendEl-mybox().HMNtrack->PositionStart[3]);
      Float_t dist1 = (utils::ArrayToTVector3(fgdTrack->PositionStart) - utils::ArrayToTVector3(mybox().HMNtrack->PositionStart)).Mag();
      Float_t dist2 = (utils::ArrayToTVector3(fgdTrack->PositionEnd)   - utils::ArrayToTVector3(mybox().HMNtrack->PositionStart)).Mag();
      Float_t prodEl1 = utils::ArrayToTVector3(fgdTrack->DirectionStart)*utils::ArrayToTVector3(mybox().HMNtrack->DirectionStart);
      
      output().FillVectorVar(FGDElDeltaTimeIniSelMu,deltainimuiniEl);  
      output().FillVectorVar(FGDElDeltaTimeEndSelMu,deltaendmuiniEl);  
      output().FillVectorVar(FGDElDistance,TMath::Min(dist1,dist2));
      output().FillVectorVar(FGDElMuonangle,prodEl1);

    }
    // Direction and position  of the segment
    output().FillMatrixVarFromArray(FGDElDir,     fgdTrack->DirectionStart,3);
    output().FillMatrixVarFromArray(FGDElPosStart,fgdTrack->PositionStart, 3);
    output().FillMatrixVarFromArray(FGDElPosEnd,  fgdTrack->PositionEnd,   3);
 
    if( track->TrueTrack ){
      output().FillVectorVar(FGDElTId,track->TrueTrack->PDG);
      output().FillVectorVar(FGDElIndex,track->TrueTrack->ID);       
      output().FillVectorVar(FGDElParentTId,track->TrueTrack->ParentPDG);
      output().FillVectorVar(FGDElGParentTId,track->TrueTrack->GParentPDG);
      if( track->TrueTrack->TrueVertex )
	output().FillVectorVar(FGDElVId,static_cast<AnaTrueVertex*>(track->TrueTrack->TrueVertex)->RooVtxIndex);
    }
    output().IncrementCounter(NFGDEl); 
  }

  // All FGD tracks. 
  AnaTrackB* allFgdTracks[NMAXTRACKS];
  int nFGD     = anaUtils::GetAllTracksUsingFGDAndNoTPC(input().GetBunch(), allFgdTracks);
  for (Int_t i=0;i<nFGD;i++){
    AnaTrackB* track = allFgdTracks[i];
    AnaFgdTrack *fgdTrack = static_cast<AnaFgdTrack*>(track->FGDSegments[0]);
    if( !fgdTrack ) continue; 
    ///checking times
    Float_t timeini = fgdTrack->PositionStart[3];
    Float_t timeend = fgdTrack->PositionEnd[3];
    output().FillVectorVar(FGDSecTimeIni,timeini);  
    output().FillVectorVar(FGDSecTimeEnd,timeend); 
 
    if( track->TrueTrack ){
      output().FillVectorVar(FGDSecTId,track->TrueTrack->PDG);
      output().FillVectorVar(FGDSecIndex,track->TrueTrack->ID); 
      output().FillVectorVar(FGDSecParentTId,track->TrueTrack->ParentPDG);
      output().FillVectorVar(FGDSecGParentTId,track->TrueTrack->GParentPDG);
    }
    if( mybox().HMNtrack ) {
      Float_t deltainimuini = TMath::Abs(timeini-mybox().HMNtrack->PositionStart[3]);
      Float_t deltaendmuini = TMath::Abs(timeend-mybox().HMNtrack->PositionStart[3]);
      output().FillVectorVar(FGDSecDeltaTimeIniSelMu,deltainimuini);  
      output().FillVectorVar(FGDSecDeltaTimeEndSelMu,deltaendmuini); 
    }
    // Direction of the segment 
    output().FillVectorVar(FGDSecCosTheta,(Float_t)(utils::ArrayToTVector3(fgdTrack->DirectionStart).CosTheta()));
    output().FillMatrixVarFromArray(FGDSecDir,     fgdTrack->DirectionStart,3);
    output().FillMatrixVarFromArray(FGDSecPosStart,fgdTrack->PositionStart, 3);
    output().FillMatrixVarFromArray(FGDSecPosEnd,  fgdTrack->PositionEnd,   3);
    output().FillVectorVar(FGDSecContained,        fgdTrack->Containment);
    output().FillVectorVar(FGDSecPiPull,           fgdTrack->Pullpi); 
    output().FillVectorVar(FGDSecAvgTime,          fgdTrack->AvgTime); 
    output().FillVectorVar(FGDSecLength,           fgdTrack->X); 
    //output().FillVectorVar(FGDSecLength,           fgdTrack->Length); 
    output().IncrementCounter(NFGDSec); 
    
  }

  // Fill ME variables. 
  for (Int_t i=0;i < mybox().nMichelElectrons;i++){
    AnaFgdTimeBin* me =  static_cast<AnaFgdTimeBin*>(mybox().MichelElectrons[i]);
    if( me->NHits[0] != 0 ) {
      output().FillVectorVar(MENHits,(Int_t)me->NHits[0]);
      output().FillVectorVar(MERawCharge,me->RawChargeSum[0]);
    }
    else {
      output().FillVectorVar(MENHits,(Int_t)me->NHits[1]); 
      output().FillVectorVar(MERawCharge,me->RawChargeSum[1]);
    }
    output().FillVectorVar(MEMinTime,me->MinTime);
    output().FillVectorVar(MEMaxTime,me->MaxTime);
    output().IncrementCounter(NME);
  }

  // Fill TPC1 track variables 
  for( Int_t i = 0; i < mybox().nTPC1tracks; i++ ) {
    AnaTrackB *track = mybox().TPC1tracks[i];
    if( !track ) continue;
    Float_t phi=TMath::ATan2(track->DirectionStart[1],track->DirectionStart[0]);
    output().FillVectorVar(TPC1TrackPhi,phi );
    if( track->TrueTrack ){
      output().FillVectorVar(TPC1TrackTId,track->TrueTrack->PDG);
      output().FillVectorVar(TPC1TrackParentTId,track->TrueTrack->ParentPDG);
      output().FillVectorVar(TPC1TrackGParentTId,track->TrueTrack->GParentPDG);
      if( track->TrueTrack->TrueVertex )	
      	output().FillVectorVar(TPC1TrackVId,static_cast<AnaTrueVertex*>(track->TrueTrack->TrueVertex)->RooVtxIndex);
    }
    output().FillMatrixVarFromArray(TPC1TrackPosStart,track->PositionStart,3);
    output().FillMatrixVarFromArray(TPC1TrackPosEnd,  track->PositionEnd,  3);
    output().FillVectorVar(TPC1TrackMom, track->Momentum);
    output().IncrementCounter(NTPC1Track);
  }

}

//********************************************************************
void numuCCMultiPiAnalysis::FillToyVarsInMicroTrees(bool addBase){
//********************************************************************

  // Variables from the numuCCAnalysis analysis (including the ones in baseAnalysis by default, otherwise addBase should be false
    if (addBase) _numuCCAnalysis->FillToyVarsInMicroTrees(addBase);

  // Fill here Variables specific for the numuCCMultiPi analysis
}

//********************************************************************
bool numuCCMultiPiAnalysis::Process(){
//********************************************************************

  //------ This function is called for each event (bunch) and each toy experiment  ----------
  return sel().GetSelection("kTrackerNumuCCMultiPi")->Apply(*_event,box()); 
}
